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A minimalist simulation model for lipid bilayers is presented. Each lipid is represented by a flexible 
chain of beads in implicit solvent. The hydrophobic effect is mimicked through an intermolecular 
pair potential localized at the "water" /hydrocarbon tail interface. This potential guarantees realistic 
interfacial tensions for lipids in a bilayer geometry. Lipids self assemble into bilayer structures 
that display fluidity and elastic properties consistent with experimental model membrane systems. 
Varying molecular flexibility allows for tuning of elastic moduli and area/molecule over a range of 
values seen in experimental systems. 



Lipid bilayer biomembranes are of fundamental impor- 
tance in cellular biology, and model membrane systems 
are fascinating physical systems in their own right. Simu- 
lation models for lipid bilayers have been developed over 
a range of resolutions (from fully atomistic descriptions 
to continuous elastic sheets) to address the many different 
length scales relevant to biological function and experi- 
mental study. The "mesoscopic" regime (~ 1 — lOOnm) 
is recognized as particularly relevant to the biophysical 
understanding of membrane systems |l(. 

Several coarse-grained bilayer models have been devel- 
oped with the mesoscopic regime in mind 0, EI 01 IE ■ 
Most of these depend upon explicit solvent to enforce 
bilayer stability. Coarse-grained solvent models do not 
provide detailed insight into the hydrophobic effect; the 
models are far too crude. Rather, such models provide 
a convenient means to enforce a bilayer stabilizing in- 
terfacial tension between solvent and lipid hydrocarbon 
tails (at considerable computational expense). Control- 
ling the interfacial tension directly, without recourse to 
explicit solvent, would seem (if possible) a more direct 
route to the same end. A few solvent-free models 0, 0, @ 
for bilayers do exist in the literature, but none include 
internal degrees of freedom for the lipids. Consequently, 
these models are unable to predict how molecular struc- 
ture influences membrane properties, phase behavior or 
realistic consequences of membrane heterogeneity. Such 
questions require flexible lipids to achieve even a quali- 
tative level of understanding. Q, 0, Q 

This letter presents a solvent-free lipid model that 
preserves the physics of lipid flexibility and hydropho- 
bic attraction. The physical properties of the studied 
membranes closely resemble those of a solvated model 
with similar lipid resolution. Q These results suggest that 
implicit solvent models may be appropriate for a wide 
class of problems in membrane biophysics. In particu- 
lar, the computational simplicity of the present model 
makes it very attractive for future studies of heteroge- 
neous bilayers, phase behavior and related phenomena 
dependent on mesoscale lipid structure. The elementary 



FIG. 1: Left: Definition of parameters used in model. Right: 
Sample conformation of tensionless membrane with 800 lipid 
molecules {cbend = 7e). Polar head beads are black, interface 
beads are gray, and hydrophobic tail beads are white. Simula- 
tions presented in this work employ 5 bead lipids exclusively. 
Modification to longer lipids with more tail beads is possible 
and straightforward. 



approach adopted for mimicking hydrophobic attraction 
holds promise for extension to lipid models beyond those 
studied here. 

Individual lipids are represented as semi-flexible chains 
of five beads (Fig. QJ. Bead 1 is identified as the polar 
head group, bead 2 is associated with the interface be- 
tween polar and hydrophobic components and beads 3-5 
comprise the hydrophobic tail region. Bonded bead-bead 
distances are constrained to have length a and bond an- 
gles are subject to a bending potential equivalent to that 
employed by Goetz and Lipowsky |3j: 



Ubend(9) = Cbend COS 6, 



(1) 



where is one of three bond angles on the molecule (Fig. 
^1 and Cbend is a positive energetic constant. There is no 
energetic cost for dihedral rotations. 

Individual beads interact through a combination of 
three different pair-potentials: 
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where c core , Ct a u and Ci n t are all positive energetic con- 
stants. With the exception of intra-molecular bead pairs 
separated by less than three bonds, the repulsive core 
interaction acts between all bead pairs and the tail dis- 
persion attraction acts between all tail-interface and tail- 
tail pairs. The soft interfacial attraction (discussed later) 
acts between all interface-interface pairs. The potentials 
are truncated at distances of 2a, 2a and 3a for the core, 
tail, and interfacial interactions respectively, and shifted 
to insure continuity of the potential. Truncation of the 
otherwise long ranged Ui n t is an essential component of 
the interaction and should not be viewed as an approxi- 
mation to a true r~ 2 potential. 

The results described below were obtained using the 
values fcsT = 0.9e, c core = 0.4e, Ct a u — 1.0c, and 
Cint = 3.0e with Cbend varied between 5.0 — 10. Oe. Low 
values of Cbend (< 5e) resulted in bilayers with a ten- 
dency to form pores and high values of Cbend (> 15c) 
gave rise to bilayers with ordered structure. The reduced 
units are calibrated by mandating that the simulations 
are conducted at 300K and that the largest observed 
area per molecule (at Cbend = 5.0c) corresponds to about 
0.7nm 2 . This results in the unit scale e = 2.75kJ/mol 
and a = 0.75nm. 

Simulations were carried out by Metropolis Monte 
Carlo using standard moves for short chains [j| and an 
additional move that attempted translation of an en- 
tire molecule. Stability was verified by bilayer assem- 
bly at constant box dimensions from a random config- 
uration of 128 molecules (Fig. |2J. All other simu- 
lations were conducted in the constant vanishing ten- 
sion/constant volume ensemble 0,0- Crystalline bilay- 
ers of 800 molecules were allowed to equilibrate prior to 
data collection, and fluidity was verified by lateral diffu- 
sion. A density of 0.07 lipids/er 3 was used throughout, 
which prevents the membrane from interacting with its 
periodic images in the z direction. During the course of 
simulation, lipids occasionally leave the bilayer to explore 
the box and later reenter the bilayer - i.e. at equilibrium 
monomers and bilayer lipids exchange ( - 3 % monomers 
depending on c be nd )■ 

In membranes with Cbend — 5 — 10c, molecular area 
scales inversely with molecular rigidity, but a simulta- 
neous increase in bilayer thickness preserves membrane 
volume (Fig. 01 ■ The thickness of a leaflet is de- 
fined by the average z distance between a molecule's 
interface bead and its final tail bead: (l z ) = ((rg — 
r%) ■ z). Unsurprisingly, stiff cr chains offer greater re- 
sistance to compression in length, resulting in thicker 
membranes (Fig. EJ. In model membrane systems, 
zero tension areas per molecule range from about 0.596 
nm 2 for DMPC(dimyristoylphosphatidylcholine]to 0.725 
nm 2 for DOPC(dioleoylphosphatidylcholine).[ll| Using 



our model, we achieve a 20% range ( 0.57nm 2 — 0.68nm 2 ) 
in areas by adjusting the chain stiffness alone. 
A linearly elastic sheet can be described by 



j = k A (L 2 -A )/A 



(5) 



where 7 is the surface tension, k A is the compressibility 
modulus, L 2 is the projected area, and Ao is the zero 
tension area. Although the simulation algorithm main- 
tains a constant thermodynamic tension 7, 7 represents 
the mechanical surface tension, which is measurable via 
the Virial stress tensor |l2j and fluctuates throughout 
the simulation (with sufficient averaging, (7) = 7). We 
measure k A by expressing Eq. [3] as |l3j 



k A = 



(7L 2 )(L 2 )-(7)(L 4 ) 
((L 2 ) 2 ) - (L 2 ) 2 : 



(6) 



and evaluating the averages over the course of our vanish- 
ing tension simulations. k A values range from 5 ± 4e/cr 2 
for Cbend = 5e to 28 ± 9e/a 2 for Cbend = 10e. Given 
our unit calibration, these values correspond to 40-224 
mJ/m 2 , in good agreement with single component phos- 
pholipid bilayers, which typically range from 60 — 270 
mJ/m 2 [|,[ii|. In contrast to the results of Refs. [l5L fl7| . 
we found that k A measurements for systems with 800 
molecules and 128 molecules agreed within error bars, al- 
though the smaller system measurements converged far 
more quickly. 

While k A provides a direct link to experiment, more de- 
tailed microscopic information is obtained by measuring 
the stress profile across the bilayer |3j, [l8| . Defining the 
local mechanical tension as a function of displacement z 
relative to the bilayer center of mass 7(2;) = P n (z) — Pt(z) 
(the difference between normal and tangential pressures), 
we measure the stress profile for systems with N=128 and 
N=800 molecules (Fig. |3|inset). The profiles agree quali- 
tatively with those obtained from fully atomic models |Tsj 
and nearly quantitatively with those obtained from sol- 
vated membranes also composed of five bead chains 0]. 
The peaks of high positive tension correspond to the po- 
sitions of the interface beads, indicating that these beads 
are holding the bilayer together, against the lateral repul- 
sions of polar heads and third and fourth beads. Undu- 
lations significantly smooth out the profile, even in mod- 
erately sized bilayers with 800 molecules. 

In our model, a strong attraction between interface 
beads mimics the hydrophobic effect. Since all "solvent" 
effects of our model are incorporated within Uint, the ef- 
fective interfacial tension is given by the interfacial con- 
tribution to the Virial tension: 
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where the sum is over all distinct pairs of interface 
beads. The factor of 1/2 corresponds to the two inter- 
faces present in a bilayer and T is thus defined in the usual 
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FIG. 2: Self-assembly of a bilayer patch of 128 lipids (fibund = 7e) in a box with constant dimensions L x — L y = 8.6<r, L z = 25a 
and periodic boundary conditions. The chosen area corresponds to that assumed by a pre-assembled bilayer at zero tension. 
Each Monte Carlo (MC) step includes (on average) an attempt to translate each bead in the simulation. 




7e 8e 



FIG. 3: Projected area per molecule (top) and leaflet thick- 
ness(bottom) as a function of molecular bending coefficient 
Cbend for membranes under zero tension. The volume (Ei z ) is 
insensitive to Cbend- Lines are to guide the eye. 



sense of the interfacial tension |20j . The surface pressure 
for each leaflet is given by the difference between total 
and interfacial tensions II = V — ( ! y)/2 so that II = T in 
the zero stress state simulated here. We measure values 
of T = 4 - 8e/cr 2 (32 - 65 mJ/m 2 ) (Figffl. Theoretical 
estimates range from 20 — 50 mJ/m 2 |7l Il9l l20j). 

The functional form chosen for Ui nt in our model is 
empirical, and was identified through trial and error mo- 
tivated by our previous experience with rigid lipid models 
0. The approximate magnitude of Ci n t given this func- 
tional form is dictated by physical necessity: Cmt = 3.0e 
leads to a stable fluid phase for a variety of Cbend values 
and physically reasonable interfacial tensions. Connec- 
tion between Ui n t and the hydrophobic effect is estab- 
lished solely on this basis. 

Membranes display obvious long wavelength fluctua- 
tions (Fig. ^) that we have quantified via the fluctuation 
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FIG. 4: Interfacial tension as defined in Eq. Error bars are 
smaller than symbol size and line is to guide the eye. Inset: 
Stress profile for systems with 128 molecules (solid) and 800 
molecules (dashed) showing the same pattern of peaks and 
valleys as observed in a similar solvated modelQ- Profiles 
correspond to systems with c bend = 7e. 



spectrum.0, @]- I n the constant zero tension ensemble, 
a flexible fluid sheet is expected to display undulations 
consistent with IfJ: 



(\Hp)\ 2 ) 



k B T (£ 4 ) 
k c (2ir\p\Y 



(8) 



the 



where k c is the bending rigidity, h(x,y) is 
local height of the membrane midplane, h(p) = 
1/L j dxdyh(x,y)e l27T P' r / L and the components of p are 
( ± 1,±2... ± n/2). The inferred rigidities are shown in 
Fig. [5] with a representative spectrum shown in the in- 
set. Wavelengths shorter than the membrane thickness 
clearly do not follow Eq. |H| and have been excluded from 
the fit as discussed extensively in prior work 0, 0, ^| . 
Bending rigidities range from k c = 2.5 — 16e (approx. 1 — 
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FIG. 5: Membrane bending rigidity (fc c ) as a function of 
molecular rigidity (cf, en d)- Data points represent a fit to the 
four longest wavelength modes; error bars represent the stan- 
dard error in the fit. Each data point corresponds to about 
one week of computation on a 2.3 GHz Athlon CPU. INSET: 
Spectrum, for Cbmd = 7e, where the line is a fit to Eq. Elwith 
k c = 13. 9e. 



8*1CP 20 J). The more flexible systems are consistent with 
experimental measurements of DGDG (digalactosyldia- 
cylglycerol), while the stiffer systems agree well with 
measurements of DLPC(dilauroylphosphatidylcholinc) 
and DMPC.[2lJ Stiffer molecules lead to stiffer mem- 
branes; this can be partially attributed to the increase in 
both compressibility modulus and membrane thickness. 
It has been suggested that one should generally expect 
k c = ItAcP/b, where d is the bilayer thickness and b is a di- 
mensionless constant QUID]. We measure b ~ 60 for the 
model presented here, well within the limits seen in other 
simulation models (b ~ 4 to b ~ 100) 0, El ES E3 
and experiment. 

Mesoscopic models provide a link between atomic level 
detail and macroscopic physical properties. The present 
description incorporates a level of realism previously lack- 
ing in implicit solvent models for lipid bilayers and should 
allow for detailed studies of biophysical questions where 
solvated models are computationally prohibitive. Addi- 
tionally, the physical picture afforded by this model is 
very much in the spirit of analytical theories that seldom 
consider water explicitly. Modeling at this level of detail 
should provide a critical link between experiment, theory 
and atomistic simulations. The simulation of inhomo- 
geneous membrane surfaces (with multiple lipid species 
and protein inclusions) is especially promising and is cur- 
rently under investigation. 
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